Proximity Effect and Multiple Andreev Reflections in Diffusive 
Superconductor-Normal-Metal-Superconductor Junctions 
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We present a theory of the current-voltage characteristics in diffusive superconductor-normal- 
metal-superconductor junctions. By solving the time-dependent Usadel equations we are able to 
describe the phase-coherent transport for arbitrary length of the normal wire. We show how the 
interplay between proximity effect and multiple Andreev reflections gives rise to a rich subgap 
structure in the conductance and how it is revealed in the non-equilibrium distribution function. 

PACS numbers: 74.45.-|-c,74.50.-|-r,73.23.-b 



Proximity effect is the modification of the properties 
of a normal metal (N) in contact with a superconductor 
(S) and it has been extensively studied in diffusive hy- 
brid nanostructures E| . Both equilibrium and trans- 
port properties 0, ^of diffusive SN systems are now 
well understood in the framework of the Usadel equa- 
tions 0. The transport through an SN interface is me- 
diated by Andreev reflection, where an electron coming 
from N with energy e below the superconducting gap 
A is converted into a reflected hole, thus transferring 
a Cooper pair to the S electrode. The time-reversed 
states involved in this process are coherent over a dis- 
tance Lc — inm{ y^M)/e, Lcj,), where D is the diffusion 
constant of N and L^f, is the phase coherence length. 

In an SNS junction the transport at finite bias is dom- 
inated by multiple Andreev reflections (MARs). Here, 
successive Andreev reflections at both S electrodes lead 
to a progressive rise of the quasiparticle energy. This pro- 
cess continues until the quasiparticle energy exceeds the 
gap energy. A microscopic theory of MARs has emerged 
recently and has been shown to describe quantitatively 
the current- voltage (TV) characteristics 0, the noise Q 
and the supercurrent Q in atomic point contacts. 

The interplay between proximity effect and MARs in 
diffusive SNS systems gives rise to a rich variety of physi- 
cal phenomena. For instance, the conductance exhibits a 
very peculiar subgap structure |^, 0| . This interplay 
is also revealed in the noise [ifj, [l^ or in the Shapiro 
steps '13]. The understanding of these experiments is a 
basic open problem in mesoscopic superconductivity. 

The theory of dissipative transport in diffusive SNS 
junctions has mainly been developed in two limits. The 
first one is the incoherent regime, when the normal metal 
length L > Lc- In this case there is no proximity effect 
and the transport can be described in terms of a semiclas- 
sical kinetic equation for the distribution function [l^ . 
This function was actually measured in Ref. ll^L and it 
was shown to exhibit a step-like structure, which is a 
manifestation of MARs. On the other hand, in short SNS 



junctions, when L <C = y^hD/A, where ^ is the su- 
perconducting coherence length, i.e. when the Thouless 
energy ex — HD/L^ exceeds A, the MARs are fully co- 
herent. In this regime the transport can be described by 
averaging the single-channel point-contact results with 
the bimodal transmission distribution for diffusive sys- 
tems 10. For the intermediate regime ^ < L < L^, 
when the interplay between proximity effect and MARs 
takes place 17j, there is no satisfactory theory so far. 

In this Letter we study the phase-coherent transport 
in diffusive voltage-biased SNS systems. We have solved 
the time-dependent Usadel equations, which allows us to 
calculate the I-V characteristics for the whole range of 
lengths from the short junction limit (L ^ ^) to the in- 
coherent regime (L ^ £,). We show that the interplay 
between proximity effect and MARs gives rise to a rich 
structure in the conductance in good agreement with ex- 
isting experiments. We also predict how this interplay is 
manifested in the quasiparticle distribution function. 

We consider an SNS junction, where N is a diffusive 
normal metal of length L < coupled to two identical 
superconducting reservoirs. We assume the SN interfaces 
to be fully transparent and neglect the suppression of the 
pair potential in the S leads near the interfaces. Our goal 
is the calculation of the current when a constant voltage 
V is applied. For this purpose we use the quasiclassi- 
cal theory of superconductivity for diffusive systems j^. 
This theory is formulated in terms of momentum aver- 
aged Green functions G(R, t, i') which depend on posi- 
tion R and two time arguments. These propagators are 
2x2 matrices in Keldysh space ("), where each entry is 
a 2 X 2 matrix in electron- hole space ('): 



G = 



G^ 



gR jrR. 
j:R gR 



(1) 



The Green functions for the left (1) and right (r) leads can 
be written as G/,^(t, t') = e-ii^i.riHlhQ^^^^^_^iyitii,^t' f-,lh ^ 

where for the chemical potentials we use /i; = 0, /i^ = eV . 
Here, Go(t) is the equilibrium bulk Green function of 
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a BCS superconductor. We transform to energy repre- 
sentation, in which the propagator G(R, e, e') depends 
on two energy arguments. It satisfies the non-stationary 
Usadel equation, which in the N region reads 



HD 



— V (G o VG) + efgG - Gfae' = 0, 



(2) 



where 73 is the Pauh matrix in electron-hole space. The 
convolution product o is defined as (A o B)(e, e') — 
J dei A(e, ei)B(ei, e'). Equation ^ is supplemented by 
the normalization condition G o G = —7T^S{e — e')I. 

In order to solve the Usadel equation it is convenient 
to use the time-dependent Riccati parametrization |l8l |. 
In this method the Green functions are parametrized in 
terms of scalar retarded (R) and advanced (A) coher- 
ence functions 7^''^(R, e, e') and 7'^'^(R, e, e'), and two 
distribution functions x(Tl, e, e') and x(R, e, e') as follows 



G^' = -2Tri N"^ o A/" o N- 
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where the inverse is defined via the o-operation. Using 
fundamental symmetries, all the Green functions can be 
obtained from 7^ and x. The transport equations for 
these functions in the N wire are 



ITT 



£ O 7^ + 7-^^ O f 



(4) 



dix - (9,7'') o — o (9,7 

ITT 

-{d,.x) o — o (9,7^) 

ITT 



[dzT ) — (a,a; 

ITT 

£ O X — X O £ 



(5) 



where £{e,e') = e-6{e — e'). Here, < z < 1 is the dimen- 
sionless coordinate which describes the position along the 
N wire. The expressions for , , and are ob- 
tained by comparing Eq. with Eq. JSJ. A solution of 
these equations can be found using the Ansatz: 



7^(e,e') 



E 



(6) 



where LOm = 2r7ieU. Other Fourier components are de- 
fined via 7^m(e) = 7o^m-n(^ + ^n)- Using the above 
Ansatz for any two functions A and _B, the Fourier com- 
ponents of {A o B){e,e') are given by [A o i3]„^„i(e) = 
An^i{e)Bi^jn{e)- The equations for the Fourier com- 
ponents "fnmi^) ^"^^ Xn,m{f) are the same as Eqs. (|4I5II . 




FIG. 1: (Color online) Normalized density of states in the 
middle of the wire as a function of energy for different wire 
lengths. Inset: minigap Ag as a function of the wire length. 



(3) where the o-product denotes now a matrix product in 
the Fourier indices and £n,m{^) = + ^n) • <^n,m- The 
boundary conditions for Eqs. (|4I5() can be easily deduced 
from the expressions of the bulk Green functions ji^l • In 
general, these equations have to be solved numerically. 

The above Ansatz leads to a current of the form I{t) = 
X^m ^™ exp(zma;ji), where wj = 2eV/fi is the Josephson 
frequency, and the current components can be written in 
terms of the Fourier components of the Green functions, 
' G„^m(e). We concentrate here on the analysis of the dc 

current /q, which we shall denote as /. It reads 



G 



N 



J^'^\eJ2Tr{r4God,GYl^ie)}, (7) 



where Gat is the normal state conductance. Next we 
express / in terms of the distribution function. It is 
possible to relate the component to the electron 

distribution function /(R, e, e') via the relation 2/ — 
) 1 - Er=o(7'" ° 7^)" ° ° (7^ ° 7^)"- Combining this 



with fundamental symmetries of the Green functions we 
can write the current as 

I = — n * E { ° ^-/]™.™ - Re [5 o /]„^ } , 

(8) 

where V = 1/2 + {g^^ o - o T^)/2tt^ and S = 
{gR o + J^^ o 9,j£"^)/7r2. Here, V describes the 

renormalization of the diffusion constant and S describes 
the spectral supercurrent. 

It is instructive to first have a look at the zero-bias 
density of states (DOS). In Fig. ^ we show the DOS in 
the middle of the wire for different values of L. The 
most prominent feature is the presence of a minigap Ag 
which scales with the wire length as shown in the inset of 
Fig. In the limit of long wires we obtain Ag « 3.2eT- 
It is worth remarking that while the DOS changes along 
the wire, the minigap is independent of the position. 
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FIG. 2: (Color online) Zero-temperature dc current as a 
function of the voltage for different wire lengths L. Upper 
inset: zero temperature critical current as a function of L. 
Lower inset: excess current as a function of L. 




. ! 




1 

L = 54 


1 . 












H — 


— — 1 — 
L = 6^ 


Ipn — 


- : : 










H — 


— — 1 — 
L = 8^ 


l] I — 
- 












H — 


— — 1 — 


1 [~1 — 










^ i 




L = 10^ 

1 





0.5 1 1.5 2 2.5 3 

eV/A 



Let us now turn to the analysis of the I-V character- 
istics. Figure |21 shows the zero-temperature I-V curves 
for different wire lengths. The main features are: (i) for 
L < ^ there is a pronounced subharmonic gap structure 
(SGS) at voltages 2A/ne (n integer). In particular, the 
curve L = 0.1^ reproduces quantitatively the result of 
the short-junction limit (ii) For L > ^ the SGS is 

progressively washed out as L increases, (iii) At high bias 
{eV >> A) there is an excess current, which is defined 
as lexc = I — GjyV . In the lower inset of Fig.|21we show 
lexc as a function of L. For L — > we recover the result 
elexc/iGN^) = 7r^/4-l of Ref.dji. In the opposite limit 
of long junctions {L ^ ^) Volkov et al. |21| found that 
lexc decays according to elexc/ {Gn^) = 0.82^/L. We 
find that lexc can be fitted to e/exc/(GivA) — 2.47^/L in 
the experimentally relevant range 3^ < L < 20^, a factor 
of three larger than in the limit of Ref. |^ 

As seen in Fig.|21 it is numerically difficult to reach the 
zero-bias limit, because the dimension of the matrices in 
Eqs. (|4I6|I scales with the voltage roughly as (2A/ey)^. 
However, the analysis of the low bias regime is not the 
most relevant. The SNS junctions usually have a negli- 
gible capacitance and their I-V curves are hence nonhys- 
teretic, exhibiting a transition from a supercurrent to a 
voltage state at the critical current. In the upper inset 
of Fig. 121 we show the value of the zero-temperature crit- 
ical current, Ic- As can be seen, the transition to the 
supercurrent branch would take place at voltages which 
are accessible to our numerical solution. 

The non-linearities in the I-Vs can more clearly be seen 
in the differential conductance, G = dl/dV, which is 
shown in Fig. |3| Notice that for L < ^ the SGS consists 
of a set of pronounced maxima at roughly eVn = 2A/n 
[n integer). For L > ^ the shape and height of the max- 
ima change drastically and new structure appears. For 
instance, the peak n = 2 appears slightly below eV = A 



FIG. 3: Zero-temperature differential conductance as a func- 
tion of the voltage for different wire lengths. The vertical lines 
indicate the position of eV = 2A/n with n = 1, 6. 



and it is accompanied by a much more pronounced max- 
imum above A. This maximum shifts towards A as L 
increases, until it merges with the peak at A. The peak 
above A is a common feature of the experimental obser- 
vations of the SGS H m Notice also that in the 
range 4^ < L < 10^ the SGS is superimposed on a back- 
ground that increases as the voltage decreases. The cor- 
rection in the conductance (as compared to Gn) dimin- 
ishes as L increases and, for instance, it reaches 15%-20% 
at low bias for L = 10^. All these features are in qualita- 
tive agreement with the experimental observations (see 
e.g. Fig. 5 in Ref. Illh . However, a quantitative compari- 
son would require us to extend our analysis to arbitrary 
transparency of the SN interfaces. 

The existence of the minigap Ag suggests that the 
SNS junction should exhibit a subgap structure similar 
to that of contacts with two different gaps 0]. Indeed, 
the conductance maximum above A appears roughly at 
eV ~ A -I- Ag. However, one must be careful with this 
analogy. First, other features expected for asymmetric 
junctions, like a peak at A — Ag, are absent. Moreover, 
our calculations show that as a function of temperature 
the conductance maxima shift following the bulk gap and, 
in particular, the peak above A survives even for tem- 
peratures higher than the minigap. The reason for the 
failure of the above argument is that at finite bias the 
minigap does not survive, and only a position-dependent 
pseudogap appears in the spectrum. 

In Fig. 21 we present our results for the distribution 
function obtained for the intermediate regime, and com- 
pare with the results for the incoherent regime of Ref. [H 
We first summarize the main results of the incoherent 
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FIG. 4: (a) dc component of the distribution function in the 
middle of the wire (z = 1/2) for three different wire lengths L 
and eV = 2.5A. The position of the chemical potentials are 
indicated by vertical lines. The dotted lines are the results 
of the incoherent model (see text), (b) dc component of the 
renormalization of the diffusion constant for the cases shown 
in panel (a), (c) The same as panel (a) for eV = 0.5A. (d) 
Spatial variation along the wire of the distribution functions 
in panel (a) for fixed energies. The curves in (a)-(d) have 
been shifted by multiples of 0.5 for convenience. 



model discussed in Ref. Il5|. The basic assumptions are: 
(i) there is no proximity effect, and (ii) the probability of 
Andreev reflection is 1 within the gap and outside. As 
a result, the dc component of / varies linearly with posi- 
tion z. This model predicts a staircase pattern, and two 
examples are shown as dotted lines in Fig. 01a,c), where 
the dc component of / in the middle of the wire is plot- 
ted. Our full calculation shows pronounced deviations 
from this staircase pattern. An example can be seen in 
the energy region — 2.25A < e < — 0.25A in Fig. ^a), 
where eV = 2.5A. In this region, the incoherent model 
predicts / = 3/4, while the full calculation gives smaller 
values in the whole window, except exactly at e = /i;. 

We explain these deviations in terms of the proximity 
effect. As shown in Fig. EJd) for eV — 2.5A, the prox- 
imity effect leads to a non-linear spatial variation of the 
distribution function /, except at the chemical potentials 
For t = \ii the distribution function varies linearly 
from 0.5 to 1, as in the incoherent model. However, for 
e = —0.75 A the proximity effect near the left electrode 
leads to an effective shortening of the wire length. This 
explains the negative deviation for the electron distribu- 
tion function in this energy range. For an energy near /i^ 
the deviation is positive, as the distribution function for 
this case (not shown here) varies from to 0.5 and the 
proximity effect takes place near z = \. 



A related aspect is the renormalization of the diffusion 
constant 2?, the dc part of which is shown in Fig.QJb). It 
reflects the enhancement of the transmission with respect 
to the normal state due to the proximity effect. Notice 
the correlation between the structure in T> and the devia- 
tions of / from the staircase pattern. In particular, at the 
chemical potentials of the superconductors 2? = 1, and / 
adopts the values predicted by the incoherent model. 

In summary, we have presented a theory of the I-V 
characteristics of diffusive SNS systems. We have stud- 
ied how the interplay between coherent MARs and the 
proximity effect is reflected in the conductance and in 
the distribution function. Our main results are: (i) we 
reproduce an additional peak above A in the conduc- 
tance, in agreement with experiments 0, llOi. llj; (ii) 
we predict the distribution function for the intermediate 
regime, which can be measured e.g. by the technique 
presented in Ref. Our work paves the way for the 
study of other transport properties in diffusive SNS sys- 
tems [Hill. 
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